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Successful prediction of ionospheric total electron content (TEC) data will 
help in correction of positioning errors in global navigation satellite systems 
(GNSS) caused by the ionosphere. This research paper proposes a prediction 
model for ionospheric TEC using a nonlinear autoregressive with exogenous 
inputs (NARX) neural network that utilizes past TEC data alongwith solar 
and geomagnetic indices namely F10.7, disturbed storm (Dst), Kp, Ap, and 
time of the day. We assess the prediction capability of our model at different 
latitudes during different solar activity years. We compare our results with 
another NARX model which uses previous TEC data along with time of the 
day, day of the year and season as exogenous parameters. The results show 
that for the solar minimum year the TEC prediction accuracy improves by 
35.71% and for the solar maximum year it improves by 31.20%. The results 
using root mean square error (RMSE), mean absolute error (MAE), 
correlation coefficient, and symmetric mean absolute percentage error 
(sMAPE) clearly indicate that solar and geomagnetic indices along with time 
of the day help in enhancing prediction accuracy of TEC across different 


latitudinal regions during both solar minimum and maximum years. 
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1. INTRODUCTION 

The signals from the global positioning system (GPS) satellites orbiting round the earth propagate 
through the ionosphere. Ionosphere is the ionized region of the earth’s upper atmosphere ranging from about 
60-1,000 km which contains significant amount of free electrons. The ionosphere is known for its dispersive 
nature and thus delays the electromagnetic signals and affects precise measurements in satellite navigation 
and communication. This influence is directly proportional to the density of electrons between the satellite 
transmitter and the receiver which quantifies the total electron content (TEC). Ionospheric TEC exhibits 
complex spatial and temporal fluctuations along with certain solar and geomagnetic related ionospheric 
disturbances [1]. The daily diurnal variations, seasonal variations, 27-day variations, 11-year solar cycle 
variations, geographic latitude, and some events like geomagnetic storms which generate sudden disturbances 
in the ionosphere are some obvious variations of ionospheric TEC [2]. Prediction of ionospheric TEC can 
thus help in improving the performance of GPS navigational systems. 

There are a number of empirical models like international reference ionosphere (IRI) [3], 
NeQuick [4], Klobuchar [5] which are available for reconstructing the state of the ionosphere. These models 
effectively predict the state of the ionosphere during quiet times but are not very accurate during periods of 
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high ionospheric activity. Several time-series models developed for prediction of TEC include fourier 
series [6], discrete cosine transform [7], and autoregressive distributed lag [8]. However, these models are not 
effective in predicting the non-stationary fluctuations present in the ionosphere. During last few decades, 
neural network techniques have been applied to predict the ionospheric TEC as well as model the various 
parameters in the ionospheric study [9]-[11]. These models consider the solar and geomagnetic factors in 
improving the prediction accuracy. Neural network models are found to outperform statistical methods. 
However, these models do not consider the temporal feature of TEC [12]. The TEC data exhibits 
non-stationarity which implies that the behavior of the ionosphere is extremely complex and nonlinear [13]. 
nonlinear autoregressive with exogenous inputs (NARX) neural network has been widely used for modeling 
non-linear dynamic systems [14]. Hence, in this paper we are proposing a NARX neural network model fed 
with time of the day, solar and geomagnetic indices as exogenous parameters to predict the TEC. Literature 
suggests number of applications modeled by NARX neural network. NARX model was tested on chaotic 
laser time series and traffic time series and it outperformed standard neural network architectures like time 
delay neural network (TDNN) and Elman architectures [14]. NARX gave better results over autoregressive 
integrated moving average (ARIMA) for wind speed prediction [15]. NARX was successfully tested for 
prediction of rainfall and showed improved results compared to linear algorithms like support vector 
machines [16]. NARX was used for traffic prediction and had superior prediction accuracy to RNN 
models [17]. The ability to predict symmetric horizontal (SYM-H) storm time index using NARX showed 
better accuracy compared to back propagation and recurrent Elman networks [18]. NARX has also been used 
for TEC prediction using time, day of the year and season as exogenous inputs [1]. 

In this study, we discuss prediction of TEC using NARX neural network model using exogenous 
parameters such as solar index F10.7 and geomagnetic indices such as disturbed storm (Dst), Kp, and Ap 
along with local time. We have studied TEC prediction at different latitudinal regions namely equatorial 
(0°-10° N), low-latitude (10°-30° N), mid-latitude (30°-60° N), and high-latitude (60°-80° N), during both the 
minimum solar activity year (2008) and maximum solar activity year (2014). We compare the results of our 
NARX model (henceforth called NARX1) with the NARX model developed for TEC prediction by 
Guoyan et al. [1] (henceforth called NARX2). 


2. METHOD 
2.1. Data and Input parameters 

The ionospheric TEC data is downloaded from JONOLAB (www.ionolab.org) [19]. To explore 
latitudinal influence on the prediction of ionospheric TEC, we selected four locations in the northern region, 
each situated at different latitudes. The selected stations are qaql (Greenland and Denmark), baie 
(Baie-Comeau and Canada), mas! (Maspalomas and Spain), and bogt (Bogota and Columbia). The latitude 
and longitude co-ordinates of these stations is shown in Table 1. 


Table 1. Selected stations in the northern hemisphere 


Region Stations Latitude and longitude 
High latitude region (60° to 80°) qaq1 (Greenland and Denmark) 60.7° N and 46.04 ° W 
Mid latitude region (30° to 60°) baie (Comeau and Canada) 49.18° N and 68.26° W 
Low latitude region (10° to 30°) mas (Maspalomas and Spain) 27.76° N and 15.63° W 
Equatorial region (0° to 10°) bogt (Columbia) 4.64° N and 74.08 ° W 


The exogenous parameters for the NARX model were selected based on the factors influencing 
TEC. We calculated the Pearson correlation coefficient for various exogenous parameters and selected the 
ones that correlate well with TEC data. Table 2 shows the values of Pearson correlation coefficient for the 
various exogenous parameters with TEC data. From the table it is clear that Ap and Kp show high positive 
correlation, Dst shows high negative correlation and F10.7 shows moderate correlation with TEC. Time of 
the day was selected since TEC exhibits diurnal characteristics. Thus, we selected solar index F10.7 and 
geomagnetic indices Dst, Kp and Ap along with time of the day. 


Table 2. Pearson correlation coefficient 
Exogenous parameters _ Pearson correlation coefficient 


Ap 0.67 
Kp 0.65 
Dst -0.63 
F10.7 0.37 
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The exogenous paramters were obtained from the space physics data facility at NASA's Goddard 
Space Flight Center, accessible through their website at https://omniweb.gsfc.nasa.gov/form/dx1.html [20]. 
The TEC data along with the solar and geomagnetic indices were sampled every one hour. We used the data 
from 1‘ July, 2008 to 26" July, 2008 (total samples=26 days*24= 624) for training the model and then tested 
the prediction capability of the model for unseen data from 27" July, 2008 to 31° July, 2008 (total samples= 
5 days*24=120). The same period was selected for training and testing the model for solar maximum year 
2014. The various exogenous parameters used are as follows: 
- Local time: TEC exhibits diurnal variation because the electron density progressively rises to a maximum 


; S P ` (2 2 
value at noon before falling to a minimum value at midnight. We compute sin (=) h and cos (=) h 


where local time, h ranges from 0-23. The sin and cos functions help in normalizing the time input in the 
range (-1,1) [1]. 

- Solar activity: solar activity significantly influences TEC values, with lower TEC during solar minimum 
years and higher TEC during solar maximum years [21]. F10.7, which correlates directly with the number 
of sunspots and solar radiation in the ultraviolet and visible spectrum, serves as a valuable predictor of 
solar activity [22]. 

- Geomagnetic activity (Dst, Kp, and Ap index): the Dst, Kp, and Ap represent geomagnetic activity. 

- Dst index: the Dst is a metric utilized to gauge the intensity of geomagnetic storms. It quantifies the 
deviation of the Earth's horizontal magnetic field component from its typical variation during calm 
days [23]. Under tranquil conditions, the Dst ranges from 0 to -50 nT, but during the occurrence of a 
powerful storm, it can surpass -200 nT. 

- Kp index: the planetary index Kp is the average standardized K-index obtained from 13 geomagnetic 
observatories between 44° and 60° in either hemisphere [21]. This parameter serves as a measure of 
geomagnetic activity. The higher the Kp index, more severe the geomagnetic activity, and greater the 
impact on the TEC. 

- Ap index: the Ap index, also known as the planetary index, is a numerical scale that ranges from 0 to 400, 
with higher values signifying more intense geomagnetic activity. This index serves as a valuable tool for 
monitoring the effects of space weather on the earth's magnetic field and assessing potential impacts on 
critical systems such as power grids, satellites, and radio communications [21]. 


2.2. Architecture of NARX1 model 

A NARX neural network is a dynamic model that incorporates recurrent feedback connections from 
the output layer to the input layer with time-delayed units. This unique architecture makes NARX networks 
an efficient tool for modeling and validation purposes, as they demonstrate better generalization and faster 
convergence compared to other conventional neural network models [24]. This can model nonlinear time- 
series systems and can be mathematically represented by: 


IE +1) = FOO, ¥E— 1)... yt -= ny) xO), x(t — 1),... x = nx)) + elt) Ep — E = + (kz +k) (1) 


where 7 (t+ 1) is the predicted TEC value at time (t +1). This depends on historical TEC values 
y(t), y(t —1),...,y(@—ny) and historical values of the independent exogenous inputs 
x(t), x(t — 1),..., x(t — nx) which in our case is time of the day, F10.7, Dst, Kp and Ap. ny and n, are the 
time delays for x(t) and y(t) respectively as shown Figure 1, e(t) is residual value at a specific point in 
time, ¢ [14]. 


Input layer Hidden layer 


x(t) 


Figure 1. NARX model with tapped delay line at the input and output, z~+ represents unit time delay 
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The NARX neural network, like a multilayer perceptron (MLP), comprises an input layer, a hidden 
layer, and an output layer in its architecture. However, the key difference is that the NARX network includes 
recurrent feedback connections from the output layer to the input layer with time-delayed units, enabling it to 
model dynamic systems and time-dependent data more effectively. The past exogenous inputs along with the 
past predicted TEC values are fed at the input of the model. The feedback connections allow the network to 
take into account the previous behavior of the system and to capture any patterns that may exist between the 
input and output variables over time. Figure 2 depicts the framework of NARX neural network. The input 
layer takes in the exogenous inputs along with past predicted TEC values. The output of the intermediate 
hidden layer H; (t) is given by: 


HAD =F |E; 9 Wp MEP) +E o Wig t-a) + a; (2) 


where, each hidden layer node is activated by the function F}. The input neuron x(¢-p) is connected to the i” 
hidden layer node with a connection weight of w,, and the output feedback neuron y(t-q) is also connected to 
the i” hidden layer node with a connection weight of Wiq. Additionally, each hidden layer node is associated 
with a bias term a;. These connection weights and bias terms are learned during the training process and they 
determine the strength and direction of the connections between the neurons in the network [25]. The 
predicted TEC value at the output } (¢+/) is given by: 


PAD =F [Ly w; Hb) + b] K 


where w;is the connection weight between the i” hidden layer node and output layer, b; is the bias of the 
network, N is the number of hidden layer nodes. Therefore, the predicted TEC value at time t+/ is a function 
of the weighted sum of the inputs passed through a nonlinear activation function and is determined by the 
weights and biases of the neural network. NARX neural network can be used in open-loop mode and closed- 
loop mode. The open-loop architecture and the closed-loop architecture of NARX neural network model is 
shown in Figures 3(a) and (b) respectively. 


Exogenous 
Inputs 


Hidden Layer Output Layer 


Feed 
Forward 


m y(t +1 
P(t + 1) Neural i 


Network 


Figure 3. NARX; (a) open-loop architecture and (b) closed-loop architecture 


In the open-loop architecture, as depicted in Figure 3(a), the NARX neural network predicts the 
output TEC } (¢+/) using historical values of both the exogenous data and the actual past TEC values. This 
configuration enables the model to utilize accurate information from the actual TEC time series, contributing 
to more precise inputs [26]. In the closed-loop configuration, as depicted in Figure 3(b), the NARX model 
feeds back the predicted TEC output to the input of the feedforward neural network. Consequently, the 
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prediction of } (¢+/) is based on the previous values of the exogenous data and the previously predicted TEC 
values. This closed-loop architecture allows the model to exploit its own predictions to make further 
estimations, potentially capturing dynamic dependencies and time-dependent patterns effectively. The 
mathematical expression for the NARX open-loop architecture is given by: 


PUD =f KO MA t-D) ny tery), MO, X-I), -Xl t-n,)) (4) 


and for the closed-loop architecture; 


PUD =H D Kt- t-ny), xD, Mt-D), ux t-n,)) (5) 


where f (.) is the mapping function of the neural network, *(¢+/) is the predicted TEC output of the NARX 
using either the past true TEC values y(t)y(t-1), wo b(tn,) or the past predicted TEC values 
I), PLD), ou D(t-ny) and the past values of the exogenous time series x(t),x(t-1)...,x(t-n,) and ny and ny are 
the input and output delays. 


2.3. Implementation of NARX1 model in MATLAB 

We have simulated the NARX model using the neural net toolbox in MATLAB. The simulated 
model is shown in Figures 4 and 5. The number of hidden layer neurons was set to 20. The delays n, and n, 
was set to two. These hyperparameters were finalized after multiple simulations. The model was run with 
different number of network feedback delays and then finalized to two where the prediction model performed 
well.The neurons in the hidden layer determine the accuracy of the model. Less neurons in the hidden layer 
fail to establish relationship of the data but increasing the number of neurons leads to more computational 
time [27]. Initially we trained the model in open loop configuration as shown Figure 4. The advantage of 
doing this is that the NARX model uses the previous values of the actual TEC time series and hence it 
receives more precise input during the training phase. Then we converted the model into closed loop for 
making predictions into the future as shown Figure 5. We perform recursive multi-step forecasting in which 
we predict one-step TEC value and use this predicted value in the next step to make next prediction. This 
process is continued till we forecast TEC values for the next 5 days. 


y(t) 
T Hidden Output 


y(t+1) 


Figure 4. NARX open loop architecture simulated in MATLAB during the training phase 


Hidden Output 


y(t+1) 


Figure 5. NARX closed loop architecture simulated in MATLAB for step ahead predictions 


Before training the network, the data was divided into three sets. 70% of the data was used for the 
training, 14% for validation, and 16% for testing. The network weights and biases were updated and the 
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gradient was calculated using the training data. The validation samples aid in fine-tuning the hyperparameters 
and the test samples help in computing the generalization performance of the model. The training algorithm 
used was Levenberg Marquardt. This algorithm updates weights and biases using less memory. Training 
stops when error starts increasing indicating that there is no further improvement in generalization. The 
output provides the predicted TEC values. The NARX model was developed using a hyperbolic tangent 
function in the hidden layer and a linear transfer function in the output layer. We calculate the root mean 
square error (RMSE), mean absolute error (MAE), correlation coefficient (r), and symmetric mean absolute 
percentage error (SMAPE) using the following formulas to evaluate the performance of the model. 


RMSE= Ezher À 


MAE= pe z fil (7) 
ie Zier Oi WiFi) _ (8) 
[eNsoi-w? ENO- 
100 GN l¥i-vi | 
sMAPE= y 9 
N “1 (yil41¥ 1/2 0) 


where f, is the predicted TEC by the model and y, is the actual TEC value, y; is the mean of the actual TEC 
values, ; is the mean of the predicted TEC values and N is the frequency of observations. 


3. RESULTS AND DISCUSSION 
3.1. Stability of the model 

After training the model we check the stability of the model by plotting the error autocorrelation 
function, error histogram, and regression coefficient. Figure 6 shows the error autocorrelation plot for qaq1 
station during the solar maximum year 2014. The error autocorrelation function is a valuable tool for 
validating the performance of a model. It helps describe the relationship between prediction errors over time, 
revealing how errors are correlated at different time lags. 


Autocorrelation of Error 1 


Correlation 


Figure 6. Error autocorrelation plot for qaq1 station during solar active year 2014 


An ideal prediction model is indicated by only one non-zero value in the autocorrelation function, 
which occurs at zero lag, signifying a perfect match between the model's predictions and the actual data [28]. 
The predicted errors are below 95% confidence limit (red dotted line) and thus are uncorrelated with each 
other. This also shows that the network has been trained satisfactorily and can be used for prediction. The 
error histogram plot in MATLAB is a graphical representation of the distribution of errors in a dataset. The 
plot shows the frequency of occurrence of errors, binned into intervals, with the x-axis representing the 
difference between predicted and actual TEC values and the y-axis representing the instances of those errors. 
Figure 7 shows a well-behaved error histogram plot which has a roughly symmetric distribution around the 
mean error value, with very few or no outliers. The plot clearly indicates that maximum samples from the 
training, validation and test dataset have a very small error equal to 0.01736. 
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The regression plot denotes an equation between the predicted TEC (output) and true (target) value 
of TEC. Figure 8 shows the regression line plot for baie station during training, validation and test phases 
during the solar maximum year 2014. Table 3 shows the slope and intercept values for the regression plot 
shown in Figure 8. The regression equation relating the predicted and the true TEC is shown on the y-axis. 
The coefficient of the equation shows proportionality between predicted TEC and the true TEC and is very 
close to 1. The constant term represents an error which should be added to the scaled target to make it as 
close as possible to the predicted output and is very close to 0 (see Figure 8). Overall, high values of R were 
achieved, specifically 99.75% for the training set, 99.04% for the validation set, 99.16% for the test data set 
and overall 99.54%. R values greater that 95% suggest excellent fit between true and predicted values [27]. 


Figures 6-8 clearly indicate that NARX network has been trained efficiently and can be effectively used for 
prediction of TEC data. 
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Figure 7. Error histogram plot for qaq1 station during solar active year 2014 
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Figure 8. Regression line plot for baie station during solar active year 2014 


Table 3. Slope and intercept values for the regression plot in Figure 8 


Training set Validation set Test set 
Slope Intercept Slope Intercept Slope Intercept 
0.99 0.16 0.96 0.58 1.00 0.14 


The validity of our model's results could not be meaningfully compared to existing literature due to 
variations in the methodologies used in each study. These differences include the use of different sets of stations, 
varying amounts of training data, and use of different solar minimum and active years. To accurately compare our 
results and validate our model, we simulated the NARX2 model using the identical architecture and parameters as 
described by Guoyan et al. [1] using time of the day, day of the year and season as exogenous inputs. 


3.2. For the solar quiet year 2008 

We check the prediction capability of our model by testing it with unseen data from 26" July to 31* 
July, 2008. The plots in Figure 9 show the observed and the predicted TEC at Figure 9(a) qaq1, Figure 9(b) 
baie, Figure 9(c) mas1, and Figure 9(d) bogt stations respectively during the solar minimum year 2008. The 
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graph illustrates how well the predictions align with the actual TEC values for each station during this 


specific time frame. 
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Figure 9. Observed and predicted TEC using NARX model for the period from 26" July 2008 to 31° July 
2008 for; (a) qaq1, (b) baie, (c) mas1, and (d) bogt stations during the solar minimum year 2008 


3.3. Solar maximum year 2014 
We also check the prediction capability of our model by testing it with unseen data from 26" July to 


31* July, 2014 for the solar maximum year 2014 (see Figure 10). The plots in Figure 10 show the observed 
and the predicted TEC at Figure 10(a) qaq1, Figure 10(b) baie, Figure 10(c) mas1, and Figure 10(d) bogt 
stations respectively during the solar maximum year 2014. Compared to solar minimum year 2008, there is a 
more deviation between the original and predicted TEC. This is specifically due to increased levels of 
ionization and increase solar activity and is consistent with the literature. 
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Figure 10. Observed and predicted TEC using NARX model for the period from 26" July 2014 to 31“ July 
2014 for; (a) qaq1, (b) baie, (c) mas1, and (d) bogt stations during the solar maximum year 2014 
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3.4. Comparison with NARX2 model 

To validate our model, we compared our results with NARX2. NARX1 model was implemented 
using time, solar index (F10.7), and geomagnetic indices (Dst, Kp, and Ap). The literature used NARX2 
model for TEC prediction using time of the day, day of the year and season as exogenous inputs. This model 
was implemented for China region for low and mid latitude regions for the solar maximum year 2011 and 
quiet year 2017. To ensure a fair comparison, we implemented NARX2 for our regions and also for the same 
solar minimum (2008) and active years (2014). Using the same model architecture (NARX), we wanted to 
see the improvement in prediction accuracy using solar and geomagnetic indices as exogenous parameters at 
all latitudinal regions and for both the solar minimum and active years. The RMSE, MAE, correlation 
coefficient and sMAPE values obtained for each of the station for both the solar minimum and active years 
using both NARX1 and NARX2 model are plotted in Figures 11(a)-(d) respectively. 
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Figure 11. Comparison of; (a) RMSE, (b) MAE, (c) correlation coefficient, and (d) sMAPE values for 
NARX1 and NARX2 model for qaq1, baie, mas1 and bogt stations during the year 2008 and 2014 


Both the RMSE and MAE metrics show that our NARX1 model can predict TEC with better 
accuracy during both the solar minimum and active years. The prediction accuracy for the solar minimum 
year 2008 shows an improvement of 35.71% (RMSE) and 42.14% (MAE) over NARX2. For the solar 
maximum year 2014, the prediction accuracy shows an improvement of 31.20% (RMSE) and 27.7% (MAE) 
over NARX2. A high correlation coefficient between true and predicted values using NARX1 model 
indicates that the predicted values are closely aligned with the true values. Further, low values of sMAPE for 
NARXI model indicate a higher level of accuracy and better performance of the forecasting model compared 
with NARX2. A lower sMAPE suggests that the forecasted values are closer to the actual values, with 
smaller percentage errors. As anticipated, the RMSE increases as latitude decreases, implying better 
prediction performance in higher and mid-latitude regions compared to lower and equatorial regions. The 
challenging aspect of modeling and predicting ionospheric TEC in low-latitude regions is attributed to the 
presence of the equatorial ionization anomaly [29]. This anomaly introduces complex and variable 
ionospheric behavior, making accurate modeling and prediction more difficult in these specific areas. The 
RMSE and MAE values for solar maximum year are higher compared to solar minimum year. This is to 
increased solar radiation, ionization, and disturbances in the ionosphere during solar maximum year and is 
consistent with the literature. The RMSE, MAE, correlation coefficient and sMAPE metrics show that 
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NARX1 model can predict TEC with better accuracy even for solar maximum year. According to our study 
findings, the NARX1 model consistently demonstrated better predictive performance for TEC values 
compared to the NARX2 model. 


4. CONCLUSION 

This study examined the use of NARX neural network for prediction of TEC data during the solar 
minimum and active years using time of the day, solar, and geomagnetic indices as exogenous parameters. The 
NARX model's ability to capture nonlinear dynamics, incorporate autoregressive terms, and include exogenous 
inputs make it a powerful tool for predicting TEC data. The error analysis based on RMSE, MAE, correlation 
coefficient and sMAPE indicate that use of solar index (F10.7) and geomagnetic indices (Dst, Kp, and Ap) as 
exogenous data helps in effective prediction of TEC. These parameters deviate from their normal values with 
changes in the ionospheric TEC. Any fluctuations in TEC due to disturbance in the ionosphere is quickly 
captured by these exogenous parameters and are effective in TEC prediction. But, NARX2 uses time of the day, 
day of the year and season as exogenous parameters. These parameters cannot capture the sudden ionospheric 
disturbances and are only useful in TEC prediction during quiet stable ionospheric conditions. Further, the TEC 
data, the solar and the geomagnetic data are all used in raw form without normalizing indicating the suitability 
of this technique for real time TEC prediction. The prediction of TEC at different latitudinal regions is also 
consistent with the literature and shows increasing RMSE with decreasing latitudes. The error autocorrelation 
function, error histogram and regression line plots also indicate the best fit of the model for TEC prediction. The 
RMSE between 0.8212 TECU to 4.0277 TECU for the solar minimum year and between 1.8952 TECU to 
6.9822 TECU for the solar maximum year indicate the suitability of this model and correct selection of the 
exogenous parameters for ionospheric TEC prediction. All the evaluation metrics clearly showed that the 
NARXI model outperformed the NARX2 model in accurately predicting TEC values. 
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